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Abstract 

A  non-destructive  inverse  method  is  developed  to  determine  internal  temperature  distribution  of  the  PEMFCs.  In  this  study,  the  attention 
is  focused  on  global  measurement  for  the  irregular  temperature  distribution  at  the  interface  between  the  carbon  plate  and  the  membrane 
electrode  assembly  (MEA)  based  on  the  measured  temperature  data  on  the  outer  surface  of  the  end  plate.  A  direct  problem  solver  capable  of 
predicting  temperature  distribution  in  the  solid  layers  of  the  PEMFC  under  various  conditions  is  incorporated  in  the  inverse  approach  to  provide 
temperature  solutions.  In  this  report,  a  concept  of  point-by-point  temperature  prediction  is  presented.  This  approach  is  particularly  suitable 
for  determining  irregular  temperature  distribution  that  is  difficult  to  handle  by  the  existing  polynomial-function  approach  [C.H.  Cheng,  M.H. 
Chang,  Predictions  of  internal  temperature  distribution  of  PEMFC  by  undestructive  inverse  method,  J.  Power  Sources,  in  press].  A  number 
of  test  cases  are  considered  in  this  study.  Some  irregular  temperature  functions  are  specified  and  regarded  as  exact  temperature  distributions 
to  predict.  Meanwhile,  the  influence  of  uncertainty  in  the  measured  temperature  data  on  the  outer  surface  is  evaluated. 

©  2004  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

In  the  past  several  years,  significant  progress  in  the  de¬ 
velopment  of  fuel  cell  technology  has  been  achieved  by 
an  increasing  number  of  experimental  [1-3]  and  theoretical 
[4-7]  studies.  The  experiments  can  help  determine  the  over¬ 
all  performance  of  the  fuel  cell  and  find  out  preferable  oper¬ 
ation  conditions.  The  theoretical  studies  help  understand  the 
physico-chemical  process  and  the  transportation  phenomena 
inside  the  fuel  cell  and  provide  detailed  information  which 
may  not  be  easily  obtainable  by  the  experiments. 

On  the  other  hand,  the  optimization  methods  are  gradu¬ 
ally  introduced  into  the  fuel  cell  design  phase.  For  example, 
a  nonlinear-constrained  optimization  procedure  to  maximize 
the  performance  of  the  cathode  with  interdigitated  air  chan¬ 
nels  in  a  PEMFC  was  presented  by  Grujicic  et  al.  [8].  In  Ref. 
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[8],  the  optimization  was  based  on  a  steady-state  single-phase 
electro-chemical  model  for  the  cathode.  In  the  study  of  Mo- 
hamed  and  Jenkins  [9],  a  genetic  algorithm  is  employed  to 
optimize  a  PEMFC  stack  design  by  searching  for  the  best 
configuration  in  terms  of  cell  number  and  the  cell  surface 
area.  Grujicic  and  Chittajallu  [10]  used  a  two-dimensional 
electrochemical  model  to  determine  the  optimal  design  of 
the  operational  and  the  geometrical  parameters  for  cathode 
of  a  fuel  cell. 

In  general,  the  electrical  energy  produced  is  accompanied 
by  an  approximately  equal  amount  of  thermal  energy  dissi¬ 
pated.  Therefore,  thermal  management  of  a  fuel  cell  is  of  great 
concerns  to  the  researchers.  In  order  to  ensure  efficient  ther¬ 
mal  management  for  the  fuel  cell,  it  is  required  to  monitor  the 
internal  temperature  distribution  of  the  PEMFC.  However, 
the  internal  temperature  of  the  fuel  cell  is  usually  not  eas¬ 
ily  measured,  especially  for  a  global  measurement.  To  have 
the  internal  temperature  information,  one  may  use  destruc¬ 
tive  methods,  in  which  a  number  of  temperature  sensors  are 
inserted  into  the  fuel  cell  to  measure  the  internal  temperature 
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Nomenclature 

C 

heat  capacity  (kJkg-1  °C_1) 

h 

heat  transfer  coefficient  (W  m-2  °C-1) 

H 

height  (m) 

J 

objective  function 

k 

thermal  conductivity  (W  m-1  °C_1) 

L 

length  of  fuel  cell  (m) 

NX,  NY, 

numbers  of  grid  points  in  x-,y~,  and 

NZ 

z-direction 

* 

q 

internal  heat  source  (W  m-3) 

t 

time 

T 

temperature  (°C) 

T 

simulated  experimental  temperature  data  (°C) 

7a 

ambient  temperature  (°C) 

w 

thickness  (m) 

x,y,z 

Cartesian  coordinates  (m) 

Greek  symbols 

p 

step  size 

Y 

conjugate  gradient  coefficient 

0 

exact  temperature  solution 

P 

density  (kgm-3) 

a 

experimental  temperature  uncertainty 

CO 

random  number  varied  —  1  and  1 

Subscripts 

C 

carbon  plate 

Cu 

copper  plate 

e 

end  plate 

ex 

exact 

g 

gasket 

Uj,  k 

grid  point  indices 

ini 

initial  guess 

Superscript 

n 

iteration  step 

directly.  Unfortunately,  the  destructive  methods  may  appre¬ 
ciably  disturb  the  original  flow  and  current  fields  inside  the 
fuel  cell,  and  may  also  be  possible  to  cause  leakage  problems 
of  the  fuel  and  oxidant  gases. 

In  order  to  resolve  these  problems,  Cheng  and  Chang  [11] 
introduced  the  concept  of  an  inverse  method  for  obtaining 
the  global  temperature  distribution  at  the  MEA/carbon  plate 
interface  in  the  PEMFC.  This  method  is  capable  of  predict¬ 
ing  internal  temperature  distribution  of  a  PEMFC  efficiently 
based  on  the  outer  surface  temperature  data  without  caus¬ 
ing  any  damage  to  the  fuel  cell.  However,  in  the  study,  the 
temperature  distribution  of  the  predicted  interface  must  be 
approximated  by  a  polynomial  function.  As  a  result,  if  pre¬ 
dicted  temperature  distribution  cannot  be  cast  into  a  form  of 
a  polynomial  function  accurately,  there  will  exists  a  remark¬ 


able  error  in  predictions.  Thus,  the  flexibility  of  the  inverse 
method  is  actually  limited  due  to  the  polynomial-function 
assumption. 

In  the  present  study,  the  existing  polynomial-function  ap¬ 
proach  is  modified  and  extended  to  the  applications  for  a  more 
flexible  form  of  temperature  distribution  at  the  MEA/carbon 
plate  interface.  The  temperature  distribution  is  predicted  by 
using  a  point-by-point  concept  which  is  not  limited  by  the 
mathematical  assumption  with  a  polynomial-function  form 
for  the  temperature  distribution.  Therefore,  even  an  irregular 
temperature  distribution  at  the  MEA/carbon  plate  interface 
can  be  predicted  accurately. 

The  validity  of  the  present  method  is  demonstrated  by 
dealing  with  two  cases  of  exact  temperature  function.  In  ad¬ 
dition,  the  influence  of  the  uncertainty  in  the  measured  tem¬ 
perature  data  on  the  outer  surface  of  the  end  plate  is  evaluated. 
Relative  performance  of  the  present  approach  is  demonstrated 
by  a  comparison  with  the  existing  method. 

Fig.  1  shows  the  schematic  of  a  typical  single-cell  PEMFC. 
The  PEMFC  shown  in  this  figure  is  equipped  with  a  polymer 
electrolyte  at  the  center.  The  polymer  electrolyte  is  sand¬ 
wiched  between  two  electrodes  and  two  gas  diffusion  layers 
to  form  a  MEA,  which  is  placed  between  two  carbon  plates 
having  machined  groves  that  provide  flow  channels  for  fuel 
and  oxidant  individually.  In  addition,  two  copper  current  col¬ 
lectors  are  attached  to  the  outer  faces  of  the  carbon  plates 
to  conduct  the  current.  In  general,  the  outer  surfaces  of  the 
copper  current  collectors  are  insulated  by  gasket  layers.  The 
single  cell  is  then  compressed  tightly  by  two  end  plates. 

While  the  cell  is  in  operation,  a  certain  amount  of  heat 
is  generated  by  the  electrochemical  reaction.  The  heat  gen¬ 
erated  must  be  conducted  toward  the  outer  surfaces  of  the 


©  Membrane  electrode  assembly  ©  Copper  current  collector 
©  Flow  channels  ©  Gasket 

©  Carbon  plate  ©  End  plate 

Fig.  1.  Schematic  of  a  single-cell  PEMFC  and  the  simulated  zone. 
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Temperature  prediction  points  at 


©  Carbon  plate 
©  Copper  current  collector 


©  Gasket 
©  End  plate 


Fig.  2.  Geometry  of  the  simulated  zone. 


end  plates  and  then  dissipated  to  the  ambient  air  at  tempera¬ 
ture  ra  by  convection.  Heat  transfer  coefficient  on  the  outer 
surfaces  is  denoted  by  h.  As  stated  earlier,  the  temperature 
distribution  at  the  MEA/carbon  plate  interface  is  predicted 
by  the  point-by-point  approach  based  on  the  measured  tem¬ 
perature  data  on  the  outer  surface  of  the  end  plate.  The  tem¬ 
perature  information  on  the  outer  surfaces  of  the  end  plates 
can  be  gathered  by  using  an  array  of  distributed  temperature 
sensors  installed  on  the  outer  surface  of  the  end  plate  or  by 
an  infrared  radiation  thermal  image  system.  Note  that  local 
temperature  distribution  at  the  MEA/carbon  plate  interface 
reflects  local  condition  of  electrochemical  reaction  to  a  cer¬ 
tain  extent.  Since  in  the  present  study  the  MEA/carbon  plate 
interface  temperature  is  determined  based  on  the  temperature 
data  on  the  outer  surface  of  the  end  plate,  the  simulated  zone 
is  ranged  from  the  MEA/carbon  plate  interface  to  the  outer 
surface  of  the  end  plate,  as  indicated  by  the  dashed  lines  in 
Fig.  E  Geometric  variables  of  the  simulated  zone  are  illus¬ 
trated  in  Fig.  2,  and  the  fixed  dimensions,  material  properties, 
and  surface  conditions  are  given  in  Table  1 . 


In  Fig.  2,  the  array  of  the  temperature  measurement  points 
on  the  outer  surface  of  the  end  plate  and  the  array  of  the  tem¬ 
perature  prediction  points  at  the  MEA/carbon  plate  interface 
are  both  shown.  Typically,  the  numbers  and  the  locations  of 
the  temperature  measurement  points,  the  temperature  predic¬ 
tion  points,  and  the  grid  points  used  in  the  numerical  analysis 
of  the  direct  problem  solver,  are  all  identical. 

2.  Optimization  method 

2.7.  Direct  problem  solver 

The  solid  materials  in  the  simulated  zone,  including  the 
carbon  plate,  copper  current  collector,  gasket,  and  end  plate  of 
the  PEMFC,  are  all  assumed  to  be  homogeneous,  isotropic 
mediums.  In  these  solid  materials,  heat  conduction  is  gov¬ 
erned  by  the  following  partial  differential  equation: 

pC—=kV2T  +  q *  (1) 

r  dt 


Table  1 


Fixed  dimensions,  material  properties,  and  surface  conditions  of  the  test  cases 


Carbon  plate 

Copper  plate 

Gasket 

End  plate 

L- 0.21  m 

L- 0.21  m 

L- 0.21  m 

L  =  0.21  m 

Wq  =  0.21  m 

Wcu  =  0.21  m 

Wg  =  0.21m 

We  =  0.21m 

Hq  =  0.003  m 

Hqu  =  0.003  m 

Hg  =  0.0003  m 

He  =  0.004  m 

£c  =  95  W  m“ 1  “C”1 

£Cu  =  380Wm-1  °C_1 

ig  =  0.17  W  m-1  °C_1 

/te  =  200Wm“1  °C_1 

Edge  convection 

Edge  convection 

Edge  convection 

Edge  and  outer  surface  convection 

h  =  10  Wm-2  °C_1 

h=  10  Wm-2  °C_1 

/i  =  10Wm_2oC_1 

h  =  10  Wm-1  °C_1 

Ta  =  25  °C 

Ta  =  25  °C 

Ta  =  25°C 

Ta  =  25  °C 
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where  q  ‘  denotes  the  internal  heat  source;  p,  C,  and  k  are  den¬ 
sity,  heat  capacity,  and  thermal  conductivity  of  the  individual 
solid  materials,  respectively;  and  T  is  the  temperature.  The 
present  approach  is  based  on  the  steady- state  thermal  behav¬ 
ior.  Thus,  the  three-dimensional  steady-state  heat  conduction 
equation  without  internal  heat  source  can  be  derived  as: 


d2T  d2T  9 2T  _ 

dx2  dy2  dz2 


The  boundary  conditions  associated  with  the  above  heat 
conduction  equation  are: 


(1)  Outer  surfaces  of  end  plates  or  edges  of  all  solid  layers: 

±k ^  =  h[T  -  ra]  (3a) 

on 

where  n  is  the  coordinate  normal  to  the  respective  sur¬ 
faces,  and  h  is  assigned  to  be  10Wm-2oC-1  for  a 
natural-convection  situation. 

(2)  Interfaces  between  two  solids: 


T\  =  T2 


(3b) 

(3c) 


where  the  indices  1  and  2  denote  the  two  successive  solid 
layers  in  contact. 


Eqs.  (2)  and  (3)  are  then  discretized  to  yield  a  set  of  simul¬ 
taneous  algebraic  equations  by  the  finite-difference  method. 
With  the  help  of  successive-over-relaxation  method  (SOR) 
[12],  the  numerical  solution  for  the  three-dimensional  tem¬ 
perature  distribution  at  NX  x  NY  x  NZ  grid  points  can  be  ob¬ 
tained.  Note  that  NZ  is  the  sum  of  NZq  ,  NZqu  ,  A fZg,  and  AZe. 


2.2.  Inverse  method 


In  the  present  study,  an  objective  function  ( J)  in  conjunc¬ 
tion  with  the  optimization  process  is  defined  in  the  following: 


words,  Taj  at  all  the  NX  x  NY  grid  points  at  the  ME  A/carbon 
plate  interface  are  regarded  as  the  major  variables  to  predict. 
Despite  of  longer  computation  time  consumed,  this  point-by¬ 
point  approach  is  capable  of  predicting  interface  temperature 
distribution  corresponding  to  any  sets  of  the  measured  tem¬ 
perature  data. 

Minimization  of  the  objective  function  makes  the  differ¬ 
ence  between  T^j  and  Te;j  vanish.  The  minimization  of 
the  objective  function  /  is  achieved  by  using  the  conjugate- 
gradient  method  described  by  Hanke  [13].  The  conjugate- 
gradient  method  is  used  to  evaluate  the  gradients  of  the  ob¬ 
jective  function  and  to  find  the  conjugate  directions  for  the 
updated  solutions  with  the  help  of  a  numerical  sensitivity 
analysis  proposed  by  Cheng  and  Wu  [14].  In  general,  in  a 
finite  number  of  iterations  the  convergence  can  be  attained. 

In  this  study,  the  inverse  method  is  further  modified  to 
become  compatible  with  the  point-by-point  approach.  For 
this  purpose,  let  7^.  .  (/=  1,  2,  . . .,  NX;  j-  1,  2,  . . .,  NY)  be 
the  nth  iterative  values  of  the  temperature  at  grid  point  (/,  j) 
at  the  interface  between  the  carbon  plate  and  MEA,  and  let 
the  first  search  direction  toward  the  minimization  of  J  be  the 
steepest  descent  direction  in  terms  of  the  gradient  functions 


(a)  Temperature  distribution  on  plane  atz  =  0.0015  m 


NX  NY 

J  =  E  E  j  -  nn2  (4) 

i= 1  7=1 

where  T^j  is  the  iterative  temperature  provided  by  the  di¬ 
rect  problem  solver  at  the  grid  point  (/,  j)  on  the  outer  sur¬ 
face  of  the  end  plate,  and  fQij  is  the  experimental  measure¬ 
ment  temperature  at  the  same  point  which  is  determined  by 
Tqij  =  Texij  +  crco,  where  Texij,  cr,  and  co  are  the  exact  tem¬ 
perature,  experimental  uncertainty,  and  a  random  number  var¬ 
ied  —  1  and  1,  respectively.  Note  that  when  o-  0,  the  exper¬ 
imental  temperature  distribution  measured  on  the  end  plate 
is  identical  to  the  exact  one.  In  addition,  the  MEA/carbon 
plate  interface  temperature  distribution  to  predict,  Tq  (. X ,  Y ), 
is  represented  by  a  matrix  of  NX  x  NY  discrete  values  of  Taj. 
The  NX  x  NY  values  of  Taj  at  the  NX  x  NY  grid  points,  that 
lead  to  minimization  of  the  objective  function  defined  in  Eq. 
(4),  are  determined  by  using  the  optimization  method.  In  other 


(b)  Temperature  distribution  on  plane  at y  =  0.105  m 


Fig.  3.  Grid  independence  of  temperature  solutions  obtained  by  direct  prob¬ 
lem  solver,  based  on  11  x  11  x  41,  21  x  21  x  41,  and  31  x  31  x  41  grids. 
The  exact  temperature  function  of  case  2  is  given  for  the  MEA/carbon  plate 
interface. 
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Exact 


0.2 


Polynomial-function  method  [11] 


a)  Case  1 


Present  point-by-point  method 


(b)  Case  2 


T 


T 


Fig.  4.  Comparison  between  polynomial-function  method  [11]  and  present  approach,  for  cases  1  and  2. 
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which  are  determined  by 


dJ 


NX  NY 


k=  1  q=  1 


i  =  1,2,  ...,NX;  j  =1,2,...,  NY 


where  the  terms  dTe^q/dTaj  are  referred  to  as  the  sensitiv¬ 
ity  coefficients.  The  task  of  the  numerical  sensitivity  analysis 
[14]  is  to  evaluate  the  sensitivity  of  the  objective  function 
/  with  respect  to  the  perturbations  of  Taj .  Therefore,  the 
terms  dT^q/dTaj  on  the  right-hand  side  of  Eq.  (5)  are  calcu¬ 
lated  by  introducing  small  perturbations  to  Taj  at  grid  point 
(i,  j )  at  the  predicted  interface.  The  perturbed  temperature 
on  the  predicted  interface  is  carried  out,  and  then  the  three- 
dimensional  temperature  solutions  are  obtained  by  the  direct 
problem  solver.  By  using  the  obtained  temperature  solutions, 
the  sensitivity  of  Te^q  at  point  (k,  q)  to  the  perturbation  in  Taj 
at  point  (/j)  is  determined  such  that  the  value  of  3T^ql3Taj 
can  be  calculated.  The  choice  of  the  magnitude  of  the  small 
perturbation  in  Taj  for  each  grid  point  at  the  MEA/carbon 
plate  interface  is  critical  for  the  numerical  sensitivity  analy¬ 
sis.  In  the  study,  the  perturbation  quantity  is  typically  chosen 
between  0.1  and  1.0. 


NX  NY 

jn-\~  1  ^  ^  ^  ^  r  r-pU  “T  1  /  rj-lU  “T  1  rjifl- )- 1  'T'ZZ  -\-  1 

J  —  Z ^  Z^  Lie/,7  ^0,1’  iCl,2»  '  •  •  ’  1Cl, NY^  • 
i=  1  7=1 


Based  on  the  conjugate-gradient  method,  the  values  of 
Taj  on  all  the  grid  points  are  updated  by 


rrn-\~  1  rrri  n  _zz 

Ci.j  =  U  i.i  -  PiJHj - 
i  =  1,2, ... ,  NX;  j  =  1,2,...,  NY 


where  each  of  the  search  directions  7zf  •  is  expressed  as  a 
linear  combination  of  the  steepest  descent  directions  and  a 
modified  vector.  That  is, 


_z?  _ 

71 Uj  ~ 


dJ  |  n  n  —  1 

dTaj  iJ  ’ 


i  =  1,2, ... ,  NX;  j  =1,2,...,  NY 


where  the  conjugate  gradient  coefficients  y" .  are  calculated 
by 


yn .  = 
YlJ 


(dJ/dTajf 


-i  2 


(dJ/dTajf-1 

i  =  1,2,...,  NX;  j  =  1,2, . 


. ,  NY 


The  step  sizes  frj  (7  =  1,2,...,  NX;  7=1,2,...,  NY)  appearing 
in  Eq.  (6)  are  to  be  determined.  In  theory,  the  values  of  fiij 
are  selected  to  minimize  the  updated  objective  function  Jn+X . 
With  the  help  of  Eqs.  (4)  and  (6),  Jn+X  can  be  deduced  as 


rj^n  1  'T’U  1  'T’U  1  \  P  ~\ 

1 C NX,  1  ’  1CNX, 2’  '  •  •  ’  1CNX,NYl  ~  ^ezjl 


NX  NY 


['T’U- (-1  /  rj-’fl 

L^ezJ  ^ci,l 


^l,l7rl,l’  ^Cl,2  2^1,2’  •  ’  TCl, NY 


i=  1  7=1 


P 1 , 7VF 7T " ? ,  .  .  .  ,  Tq  nx,  t  /3nx,  1 nNX,  i ,  Tq  nx  2  P NX, 2 ^NX, 2’  ■■■  ’  ^CNX,NY  &NX, NY ^NX, NY )  ^ez ,  j ] 


(9) 


By  introducing  a  first-order  Taylor  series  approximation  into 
the  above  expression,  the  (n+ 1  )th  objective  function  becomes 


NX  NY 


Jn+1  = 


i= 1  7=1  L 


ar2  .  dTn  .  ar*  . 

(*3,7  -  W -  ( +  A.2<2^  +  ■  ■  ■  +  + 


DTn 

01Cl,\ 


0Tn 

01Cl,2 


0Tn 

01C\,NY 


dTn .  ar2 .  ar2 . 

+&VX, 1  dTn&lJ  +  PNX, 2KnNXj  H - +  PnX,NY^nx,NY  ^ 


-i  2 


CM,1 


0Tn 

01 C  NX,  2 


0Tn 

U1CNX,NY 


(10) 


Making  the  derivative  of  Jn+l  with  respect  to  fiij  vanish  gives 


NX  NY 
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Eq.  (11)  represents  a  set  of  NX  x  NY  linear  algebraic  equa¬ 
tions  which  could  be  solved  simultaneously  to  yield  the  op¬ 
timal  step  size  foj  (i=  1,2,...,  NX: ; 7=1,2,...,  NY)  for  each 
TCij  at  point  (ij).  The  solution  can  be  evaluated  numerically 
by  means  of  the  Gaussian  elimination  method  or  other  similar 
solvers. 

Based  on  the  above  method,  the  temperature  distribution 
at  the  ME  A/carbon  plate  interface,  Taj,  is  updated  in  it¬ 
eration  until  an  optimal  temperature  distribution  satisfying 
the  /-minimization  criterion  is  obtained.  Typically,  the  /- 
minimization  criterion  of  the  temperature  prediction  is  set 
with  /<  1.0  x  10-5. 

Two  cases  of  known  temperature  distribution  functions 
at  the  MEA/carbon  plate  interface  are  specified  as  the  ex¬ 
act  temperature  solutions.  For  comparison,  both  the  present 
approach  and  the  existing  approach  [11]  are  applied  to 
evaluate  the  relative  performance  of  the  present  method. 
The  two  exact  temperature  functions,  0  (x,  y),  are  given 
as: 


Case  1  :  0(x,  y)  =  80  +  lOx  —  80x2  +  20y  —  80y2  —  5xy 


Case  2  :  0(x,  y) 


(12) 


(13) 


where  0  is  in  °C  and  x  and  y  are  in  m. 

Firstly,  numerical  checks  have  been  performed  to  ensure 
grid-independence  of  the  direct  problem  solver.  Numerical 
predictions  of  the  three-dimensional  temperature  distribu¬ 
tion  in  the  PEMFC  are  obtained  based  on  three  different  grid 
systems,  and  the  obtained  results  are  compared.  Three  grid 
systems,  having  NX  x  NY  xNZ=  11x11  x  41,  21  x  21  x  41, 
and  31  x31  x  41  grids,  respectively,  are  tested.  For  the  ther¬ 
mal  condition  in  which  the  MEA/carbon  plate  interface  tem¬ 
perature  is  specified  with  the  exact  temperature  function  of 
case  2,  the  numerical  temperature  solutions  on  the  planes 
at  z  —  0.0015  m  and  at  y  —  0.105  m  are  shown  in  Fig.  3.  In 
the  figure,  it  is  found  that  an  increase  in  grid  number  from 
21  x  21  x  41  to  31  x  31  x  41  produces  no  appreciable  dif¬ 
ferences  between  the  two  sets  of  solutions.  Therefore,  the 
grid  system  of  21x21x41  grids  is  used  in  this  study 
typically. 

The  two  exact  temperature  functions  defined  in  Eqs.  (12) 
and  (13)  are  introduced  individually  for  the  thermal  bound¬ 
ary  condition  at  the  MEA/carbon  plate  interface  to  obtain 
the  three-dimensional  temperature  solutions  in  the  PEMFC 
by  the  direct  problem  solver.  From  the  obtained  numerical 
three-dimensional  temperature  solutions,  the  temperature  so¬ 
lutions  on  the  outer  surface  of  the  end  plate  are  recorded  and 
regarded  as  the  exact  temperature  distribution  on  the  end- 
plate  surface  (/ex/j)-  The  exact  temperature  TQX  ij  is  added 
by  a co  to  provide  the  simulated  experimental  data  T^j  which 
are  required  for  testing  the  inverse  approach,  as  indicated 
earlier.  The  inverse  approach  is  acceptable  only  if  the  pre¬ 


diction  of  the  MEA/carbon  plate  interface  temperature  dis¬ 
tribution  based  on  the  simulated  experimental  temperature 
data  is  identical  to  the  exact  temperature  function  6  (x,  y) 
specified. 


3.  Results  and  discussion 

Fig.  4  shows  the  advantages  of  the  present  approach.  In 
this  figure,  the  predicted  MEA/carbon  plate  interface  temper¬ 
ature  distributions  by  the  present  point-by-point  method  and 
by  the  polynomial-function  method  [11]  are  compared  with 
the  exact  temperature  functions.  Plotted  in  the  left  port  are  the 
results  for  case  1,  and  the  results  for  case  2  are  plotted  in  the 
right  portion.  It  can  be  observed  that  for  case  1  of  which  the 
exact  temperature  function  to  predict  is  actually  a  polynomial 
function  expressed  by  Eq.  (12),  both  the  polynomial-function 
and  the  present  methods  lead  to  the  temperature  solutions 
identical  to  the  exact  temperature  function.  However,  for  the 
irregular  exact  temperature  function  of  case  2  given  in  Eq. 
(13),  the  polynomial-function  method  fails  to  provide  accu- 


x  or  y  [m] 

Present  point-by-point  method 

Fig.  5.  Detailed  temperature  data  along  x  -  0. 105  m  and  y  =  0. 105  m:  com¬ 
parison  between  polynomial-function  method  [11]  and  present  approach,  for 
case  2. 
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NXxNY=  21x21 


(a)  a  =  0.001  [°C] 


(b)  a  =  0.01  [°C] 


Fig.  6.  Effects  of  experimental  temperature  uncertainty  and  number  of  temperature  prediction  points  on  predictions  of  temperature  distribution  at  ME  A/carbon 
plate  interface. 


rate  predictions.  Nevertheless,  on  the  other  hand,  the  present 
approach  still  leads  to  a  satisfactory  temperature  distribution 
which  closely  agrees  with  the  exact  temperature  function. 

In  order  to  have  a  closer  look  at  the  discrepancy  be¬ 
tween  the  two  sets  of  results  obtained  by  the  two  approaches, 
the  predicted  temperature  profiles  along  the  central  line  at 
v  =  0. 105  m  and  the  central  line  at  y  =  0. 105  m  are  plotted  and 
compared  with  the  exact  profiles  for  case  2  (Fig.  5).  The  re¬ 
sults  obtained  by  the  polynomial-function  method  are  shown 
in  the  upper  plot,  and  the  results  by  the  present  approach 
in  the  lower  plot.  The  dashed  and  dash-dot  curves  indicates 
the  predicted  solutions  along  v  =  0. 105  m  and  y  =  0. 105  m,  re¬ 
spectively.  The  corresponding  exact  profiles  are  indicated  by 
the  solid  and  the  thicker-solid  curves.  It  is  clearly  observed 


that  for  this  case  the  polynomial-function  method  produces 
significant  errors  in  the  predictions,  while  the  present  ap¬ 
proach  leads  to  satisfaction. 

It  is  suspected  that  the  uncertainty  (cr)  of  the  simulated 
experimental  temperature  data  (T^j)  may  play  an  important 
role  in  the  accuracy  of  the  solution.  Meanwhile,  the  number 
of  the  temperature  prediction  points  (NX  x  NY)  on  the  end 
plate  surface  may  also  be  an  influential  factor.  Fig.  6  shows 
the  effects  of  uncertainty  and  number  of  prediction  points  on 
predictions  of  temperature  distribution  at  MEA/carbon  plate 
interface.  It  is  found  that  the  error  of  predictions  is  sensi¬ 
tive  to  the  measurement  uncertainty,  and  the  relation  between 
the  prediction  error  and  the  measurement  uncertainty  is  de¬ 
pendent  on  the  number  of  prediction  points.  When  11x11 
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Final  solution  (Exact  solution) 

(a)  Tm  =  80  °C  (constant)  (b)  TM  =  80  -  (80/0.2 1 )  x  °C 


Fig.  7.  Effects  of  initial  guess  on  the  final  solution,  for  case  2. 
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Fig.  8.  Iterating  objective  functions  with  different  initial  guesses,  for  case  2. 

prediction  points  are  used,  an  uncertainty  of  0.001  °C  in  sim¬ 
ulated  experimental  data  causes  no  appreciable  prediction 
error.  As  a  is  elevated  to  be  0.01  °C,  the  average  fluctuation 
of  the  temperature  predictions  for  the  MEA/carbon  plate  in- 


Fig.  9.  Variation  in  temperature  profile  along  the  central  line  at  y  =  0. 147  m, 
for  initial  guess  of  7jni  =  80  °C. 


terface  is  only  increased  to  approximately  0.24  °C.  However, 
when  21x21  prediction  points  are  used,  the  fluctuation  of  the 
predictions  becomes  severer.  The  average  value  of  the  fluctu¬ 
ation  reaches  2.26°Cata  =  0.01°C,  which  is  not  acceptable 
in  practical  applications.  The  increase  in  prediction  error  with 
the  number  of  prediction  points  is  expectable  in  the  present 
approach  since  in  this  point-by-point  approach,  temperature 
at  all  the  NX  x  NY  prediction  points  on  the  MEA/carbon  plate 
interface  (Taj)  are  regarded  as  the  major  variables  to  predict. 
A  larger  number  of  the  prediction  points  means  a  larger  group 
of  the  major  variables,  and  therefore,  results  in  difficulties  in 
reducing  the  magnitudes  of  the  prediction  errors.  To  improve 
the  performance  of  the  present  approach,  the  fluctuation  of 
the  temperature  predictions  must  be  further  reduced  in  future 
works. 

One  may  have  reasons  to  suspect  that  the  point-by-point 
method  might  not  lead  to  unique  solution  when  different 
initial  guesses  are  used.  To  test  the  uniqueness  of  the  pre¬ 
dictions,  two  kinds  of  initial  guess  are  adopted.  One  is  a 
uniform  temperature  distribution,  7jni  =  80  °C,  and  the  other 
is  a  linearly-varied  distribution  from  0  to  80°  C  given  with 
Tini  =  80  —  (80/0.21) x°C.  Fig.  7  shows  the  predicted  tem¬ 
perature  distributions  yielded  from  the  two  different  initial 
guesses  by  the  present  approach  for  case  2.  It  is  interesting  to 
find  that  for  this  case  only  one  unique  solution  is  obtained  re¬ 
gardless  of  the  initial  guesses,  and  both  the  obtained  solutions 
are  in  good  agreement  with  the  exact  function.  The  associated 
objective  functions  varying  with  iteration  step  are  shown  in 
Fig.  8.  The  convergence  criterion  is  set  with  J<  1  x  10“5.  It 
is  observed  that  the  objective  functions  are  decreased  rather 
rapidly.  The  objective  function  of  the  test  case  with  the  initial 
guess  of  linearly  varied  temperature  distribution  is  generally 
greater  than  that  the  uniform  temperature  distribution;  how¬ 
ever,  in  five  to  six  iteration  steps,  both  cases  achieve  conver¬ 
gence  and  lead  to  identical  solutions  by  using  the  conjugate- 
gradient  method.  For  the  test  case  with  the  initial  guess  of 
uniform  temperature  distribution,  the  variation  in  tempera¬ 
ture  profile  along  the  central  line  at  y  =  0.147  m  is  illustrated 
in  Fig.  9.  It  is  observed  that  the  temperature  profile  is  changed 
immediately  to  a  profile  which  is  close  to  the  final  solution 
at  just  the  second  iteration  step.  Then,  in  a  few  iterations  the 
solution  is  rapidly  improved  to  the  one  satisfying  the  conver¬ 
gence  criterion. 

4.  Concluding  remarks 

A  non-destructive  inverse  method  is  developed  to  de¬ 
termine  internal  global  temperature  distribution  at  the 
MEA/carbon  plate  interface  in  PEMFCs  based  on  the  mea¬ 
sured  temperature  data  on  the  outer  surface  of  the  end  plate. 
A  concept  of  point-by-point  temperature  prediction  is  pre¬ 
sented  and  demonstrated  by  a  number  of  test  cases.  Some 
irregular  temperature  functions  are  specified  and  regarded  as 
exact  temperature  distributions  to  predict.  Meanwhile,  the  in¬ 
fluence  of  uncertainty  in  the  measured  temperature  data  on 
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the  outer  surface  is  evaluated.  Relative  performance  of  the 
present  approach  is  demonstrated  by  a  comparison  with  the 
polynomial-function  method  proposed  by  Cheng  and  Chang 
[1 1]  in  predicting  the  exact  temperature  distributions. 

It  is  observed  that  for  case  1  of  which  the  exact  tempera¬ 
ture  function  to  predict  is  actually  a  polynomial  function,  both 
the  polynomial-function  and  the  present  methods  lead  to  the 
temperature  solutions  identical  to  the  exact  temperature  func¬ 
tion.  However,  for  the  irregular  exact  temperature  function  of 
case  2,  the  polynomial-function  method  produces  significant 
errors  in  the  predictions,  while  the  present  approach  leads  to 
satisfaction. 

The  uncertainty  (<r)  of  the  simulated  experimental  temper¬ 
ature  data  (TCij)  and  the  number  of  the  temperature  prediction 
points  (NX  x  NY)  on  the  end  plate  surface  are  influential  fac¬ 
tor  in  the  accuracy  of  the  solution.  It  is  found  that  the  error  of 
predictions  is  sensitive  to  the  measurement  uncertainty,  and 
the  relation  between  the  prediction  error  and  the  measure¬ 
ment  uncertainty  is  dependent  on  the  number  of  prediction 
points.  The  prediction  error  increases  with  the  measurement 
uncertainty  or  the  number  of  prediction  points. 

Effects  of  the  initial  guess  on  the  uniqueness  of  the  pre¬ 
dicted  solution  are  investigated.  It  is  found  that  for  case  2 
only  one  unique  solution  is  obtained  regardless  of  the  ini¬ 
tial  guesses,  and  the  obtained  solutions  from  different  initial 
guesses  are  in  good  agreement  with  the  exact  function. 
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